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We study the role of quantum fluctuations of atomic nuclei in the real-time dynamics of non- 
equilibrium macro-molecular transitions. To this goal we introduce an extension of the Dominant 
Reaction Pathways (DRP) formalism, in which the quantum corrections to the classical overdamped 
Langevin dynamics are rigorously taken into account to order Ti 2 . We first illustrate our approach 
in simple cases, and compare with the results of the instanton theory. Then we apply our method 
to study the C7 eq —> C7 ax transition of alanine dipeptide. We find that the inclusion of quantum 
fluctuations can significantly modify the reaction mechanism for peptides. For example, the energy 
difference which is overcome along the most probable pathway is reduced by as much as 50%. 
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I. INTRODUCTION 

Classical or ab-initio molecular dynamics (MD) sim- 
ulations have become a standard tool to investigate a 
wide range of physical systems, with widespread applica- 
tions in chemistry, material science and molecular biol- 
ogy. Such approaches are based on the assumption that 
the atomic nuclei can be treated as classical particles pQ . 

A classical description can be considered reliable for 
most atomic species comprising organic molecules and 
materials. In fact, quantum fluctuations of carbon, oxy- 
gen, nitrogen atoms at room temperature are expected to 
lead to small corrections. On the other hand, quantum 
effects are expected to play a much more important role 
in the dynamics of the lightest atomic species. For ex- 
ample, typical quantum fluctuations of a hydrogen atom 
around its equilibrium configuration in a macro-molecule 
at room temperature can be shown to be of the order of 
fractions of the Bohr radius. 

An efficient method was developed [2 to account for 
quantum effects in the evaluation of thermal averages, 
in the semi-classical and non-degenerate temperature 
regime 



p — « a 2 [fi = l/k B T), (1) 
m 

where a is a typical length-scale characterizing the in- 
teraction between atoms. In such an approach, av- 
erages of arbitrary configuration-dependent observables 
can be evaluated to order ft 2 accuracy by simply replac- 
ing the potential energy U (X) in the Boltzmann's weight 
P(X) oc e _/3C/ ( x ) with an effective semi-classical poten- 
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tial E/q(X), which reads 
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In such an expression, X = #2, . . . , xn) is a point in 
the 37V-dimensional configuration space defined by the 
set of all atomic coordinates and 



are characteristic parameters which set the length scale of 
quantum fluctuations of the particles. For example, for a 
hydrogen atom at room temperature, A = 1.3 x 10 -4 nm 2 , 
therefore y/X is about 20% of the Bohr radius. 

Clearly, this approach is only useful to investigate ther- 
modynamical properties of molecular systems. Account- 
ing for quantum corrections to dynamics and kinetics is 
in general a much more challenging task. To this goal, a 
number of methods have been proposed in the literature, 
such as centroid molecular dynamics [7 , or instanton- 
based approaches [5HT2] . 

All of these methods represent useful tools to target 
specific questions. For example, the centroid method can 
be used to investigate the real-time evolution of quan- 
tum many-body systems over short time intervals. On 
the other hand, it becomes very inefficient to investigate 
the long-time dynamics of thermally activated reactions. 
The reason is that, like any algorithm based on the inte- 
gration of the equation of motion, it wastes most of the 
computational time to simulate the exploration of the 
meta-stable states, i.e. when the system is not undergo- 
ing the transition. 

Instanton-based methods provide an elegant and pow- 
erful tool to compute quantum corrections to the reac- 
tion rates. On the other hand, they do not yield direct 
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information about the real-time non- equilibrium dynam- 
ics, since they are based on a path integral representation 
of the quantum partition function. 

In this paper, we introduce a formalism which comple- 
ments the existing methods and allows to efficiently and 
rigorously investigate the real-time evolution of (macro)- 
molecules in non-equilibrium conditions. Such a fully mi- 
croscopic approach is based on the path integral represen- 
tation of the solution of a Fokker-Planck equation which 
includes order ft 2 quantum corrections. In particular, 
our method is useful to investigate reaction mechanisms, 
since it yields a natural and unbiased reaction coordinate, 
and it allows to predict the time evolution of arbitrary 
observables during the most probable reaction pathways. 
In many cases of interest, such information can be di- 
rectly compared against experimental data. For exam- 
ple, in the context of protein folding, information about 
the reaction mechanism are available from the so-called 
phi-value analysis [14 or from single-molecule kinetic ex- 
periments — see e.g. [15] and references therein — . 

The semi-classical extension of the DRP method we 
present in this work is computationally very efficient: on 
the one hand, it avoids wasting computational time to 
simulate the dynamics when the system is not undergoing 
a transition to the final state. On the other hand, the 
inclusion of quantum corrections to order h 2 does not 
involve a significant increase of the computational cost 
of the calculation, since it requires to compute quantities 
which are already evaluated in the classical approach. 

We first illustrate our approach on a very simple two- 
dimensional toy-system. Then, we compare the effects 
of quantum fluctuations on the reaction pathways for H2 
dissociation on the Cu(llO) surface obtained in our ap- 
proach — which holds in non-equilibrium conditions — 
and in the instanton method — which applies to equilib- 
rium conditions — . 

We then perform an application to a realistic molecular 
transition: the C7 eq — >• C7 ax re- arrangement of the ala- 
nine dipeptide. This is a representative example of a bio- 
molecular conformational reaction, where the formation 
and breaking of hydrogen bonds is an important driving 
force. While quantum fluctuations of hydrogen atoms 
may play an important role in the hydrogen-bonding dy- 
namics, they are usually neglected in standard biochemi- 
cal simulations. However, as we shall see, their inclusion 
has significant effects on the reaction mechanism. 

The paper is organized as follows. In section [II] we 
discuss the real-time quantum diffusive dynamics of a 
molecular system, in the strong friction limit. Section 
ITT] represents the core of the paper, where we introduce 
the leading quantum corrections to the DRP approach. 
In section |IV[ we illustrate our method in a simple toy 
model and study the H2 dissociation in Cu. In section [V] 
we apply it to the alanine-dipeptide transition. Conclu- 



sions and perspectives are summarized in VI The math- 
ematical details of the derivations are reported in two 
appendixes. 



II. QUANTUM DIFFUSIVE DYNAMICS IN 
THE HIGH FRICTION LIMIT 

Our goal is to include quantum corrections to the dy- 
namics of systems which, in the classical limit, can be 
described by the overdamped Langevin dynamics. To 
this end, we consider the theory of quantum dissipative 
systems in the high-friction and semi-classical regime. 
Namely, if 7 is a friction coefficient which describes the 
strength of the coupling of the system to the heat-bath, 
and ro is the shortest time scale of the dynamics we are 
interested in, we consider the long-time evolution of the 
system in the limit: 



£>r > I/7, 
(3<yh <C 1. 



(4) 
(5) 



The first inequality defines the overdamped regime. The 
second condition implies that quantum coherence effects 
are negligible, see e.g. [3 and references therein. As 
a consequence, the real-time quantum dynamics in this 
limit is completely specified by the diagonal part of a 
reduced density matrix, in which the heat-bath degrees 
of freedom are traced out: 



P(X,t) = (X|Tr y p(t)|X). 



(6) 



In this definition, p(t) is the time-dependent density oper- 
ator of an enlarged Hamiltonian system which comprises 
the molecule degrees of freedom X and the heat-bath de- 
grees of freedom Y. The trace Try is performed over the 
heat-bath variables Y only. 

The equation which determines the probability density 
P(X, t) can be derived using the path integral represen- 
tation of the time-dependent density matrix. It has been 
recently shown that, in the case of a one-dimensional 
quantum particle interacting with an external potential 
U (x), the probability density P(x, t) to leading-order in A 
obeys the Quantum Smoluchowski Equation (QSE) [4] [5] 
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where D = 1/(7717/?) is the classical diffusion coefficient. 
The same result was obtained also by Coffey and cowork- 
ers, in an approach based on the thermal Wigner function 
[6]. Notice that at finite temperature, and high- friction 
regime, the leading quantum corrections are independent 
of the friction coefficient [13]. 

In appendix [A] we derive the multi-dimensional gener- 
alization of such an equation to a system of N atoms in 
contact with a heat-bath, and obtain 
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Di = 1/(771^/37) is the classical diffusion constant of the 
i-th atom of mass rm, while the functions 
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Lx(X) = (3^2 Xi V-C/(X) 



L 2 (X) = |V[/(X)| 2 
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(9) 
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account for quantum corrections and appear also in the 
leading quantum correction to Boltzmann's weight — cfr. 
Eq. ([2| — . Indeed, it is immediate to verify that 



P eq (X.) = const, e 



-pu Q (x.) 



(11) 



is the stationary solution of Eq. ([8]). 

In appendix [B] we show that the non-equilibrium prob- 
ability density P(X, t) which solves Eq. ([8| can in prin- 
ciple be sampled by integrating an associated quantum 
Langevin Eq. (QLE) with a multiplicative noise [4]: 



+v/2A(i + ii(X)) ii(t), 



(12) 



where are delta-correlated white Gaussian noises, 

with unit variance and Qi(x) are effective "quantum 
forces", whose definition depends on the choice of the 
stochastic calculus. In particular, in the so-called Ito 
calculus Qi reads 



1 - 



Qf °(X) = ^V,L 2 (X) - Li(X) Vil7(X). 



(13) 



In order to investigate the effect of the quantum cor- 
rections to the Langevin dynamics, it is instructive to 
consider the diffusion close to potential energy extrema, 
in local harmonic approximation. In this limit, the quan- 
tum Langevin Eq. ( 12 ) in the Ito calculus reads 



X 



1 



t + A/3((Tr7i )i - Ho)\ Ho (X - X ) 



t + Xf3TiH 



(14) 



where Ho is the Hessian matrix at the extremum con- 
figuration Xo. For sake of simplicity we have assumed 
that each degree of freedom is characterized by the same 
quantum parameter A. 

From Eq. (14) it follows that, near the extrema of 



the potential energy surface, the quantum contribution 
to the diffusion coefficient which multiplies the random 
force rj(t) can be re- absorbed by a rescaling of the thermal 
energy. 
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1 
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(l + A/3Tr7i ). 



(15) 
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FIG. 1: Interpretation of the quantum corrections to the 
stochastic diffusion: in the stable wells quantum fluctuations 
lead to an increase of the diffusion constant, hence the quan- 
tum diffusive motion is qualitatively similar to a classical one 
at a higher temperature. However, in the vicinity of a saddle 
where TrHo < 0, quantum effects reduce the diffusion con- 
stant. Hence, the quantum diffusive motion is qualitatively 
analog to a classical one, at a lower temperature. 



In particular, close to a minimum of the potential en- 
ergy surface one has TrHo — const. > 0, hence the quan- 
tum system diffuses like a classical one in which the heat- 
bath has a higher temperature (see left panel of Fig. [I]). 
On the other hand, near the saddles of ?7(X), i.e. in 
the transition regions, the Hessian matrix is not positive 
definite and its trace can become negative. In this case, 
the quantum system diffuses like a classical one in which 
the heat-bath has a lower temperature (see right panel 
of Fig. [Tf. 

We note that, in the particular case of one-dimensional 
systems, the quantum force Qj to (X) in local harmonic 
approximation vanishes identically. Hence, near the ex- 
trema of the potential energy, the entire o(A) correction 
to the one-dimensional Langevin dynamics is equivalent 
to a rescaling of the temperature. In higher dimensional 
systems this is in general no longer the case, since the 
o(A) correction to the force is not identically null. 

We also note that the effective lowering of the tem- 
perature induced by the quantum effects may hardly af- 
fect the analysis of thermodynamical quantities, since the 
transition regions give in general small contributions to 
equilibrium averages. On the other hand, it may have 
an important effects on non-equilibrium reactive trajec- 
tories, which by definition cross the transition region. 

In practice, for typical molecular systems, the direct 
integration of the QLE ( 12 ) can only be used for investi- 



gating very fast processes, or small thermal fluctuations 
around the local equilibrium configurations. On the other 
hand, for most molecular systems, integrating such an 
equation of motion to investigate the long-time dynamics 
of a rare activated transition would be computationally 
extremely expensive. In the next section we discuss how 
this difficulty can be rigorously overcome in the DRP 
approach. 
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III. QUANTUM CORRECTIONS TO THE 
DOMINANT REACTION PATHWAYS 

The DRP approach was originally developed to study 
the dynamics of rare thermally activated transitions in 
systems obeying the classical Langevin equation [T6HT9] . 
A remarkable advantage of the DRP approach with re- 
spect to the MD algorithm is that the computational cost 
of determining the most probable pathways in a rare ther- 
mally activated transition depends neither on the height 
of the free energy barrier which must be overcome, nor 
on the existence of gaps in the time scales associated with 
the system's dynamics. 

The DRP method has been tested so far on confor- 
mational reactions of toy-models [17] and biomolecular 
systems [20-22 . The same approach was then applied 
to investigate ab-initio both chemical reactions [23] and 
the folding of a peptide chain [24 . In these two simula- 
tions, the molecular energy U (X) and its first and second 
derivatives were obtained directly from the calculation of 
the ground state electronic structure, without resorting 
to empirical force fields. 

We now extend the DRP formalism to account for 
quantum corrections in the Langevin dynamics of the 
atomic nuclei. In appendix [B] we derive the path integral 
representation of the solution of the QSE 

P(X f ,t\Xi) = jV(X/,Xi) / 'pX e -(^//[x]+s e Q // [x]) < 

(16) 

the factor A/"(X/,Xi) is defined in the appendix [B] and 
does not affect the relative statistical weight of the reac- 
tion paths. The functional 5 e //[X] is called the (classi- 
cal) effective action and is given by 



S eff[ 

with 



y e// (x) = J2 



;V A/? 2 ( 2 



1 = 1 



Vi£/(X)| 2 - ^V^(X) ) . (18) 



Quantum effects are taken into account through the 
term 



Sf ff [X] = fdrV^lMr)}, 
Jo 



(19) 



where 



N 



K//(x) = E^ 2 i^ c/ ( x )i 2i >( x ) 

+i/?AVi • Qi(X) + ^AViLi(X) • V;t/(X). (20) 

V e ff(X.) and V®f(X) will be referred to as the classical 
and the quantum component of the effective potential, 



respectively. Note that they depend on the molecular 
energy J7(X), on the viscosity, and on the temperature 
of the heat-bath. 

For large systems, evaluating the quantum part of the 
effective potential may be quite computationally expen- 
sive, since this term contains a summation over derivates 
of the potential energy up to fourth order. A reduction of 
the computational cost can be obtained by restricting the 
summation in the quantum terms L 1 (X.) and L 2 (X) to 
the hydrogen atoms only. This is a good approximation, 
since the quantum constants of the heavier atoms in 
biomolecules is at least one order of magnitude smaller. 
In addition, if the reaction under investigation is ther- 
mally activated, the potential energy barriers which must 
be overcome are much larger than the average thermal 
energy 1//3. Thus in this case, solely the leading terms 
in the expansion of V®j(x) in powers of 1//3 may be 
retained: 



N 



V 



eff 



(X) ~ 5> 2 A 



i|V;£/(X)| 2 MX)- 



> (21) 



where the dots denote the sub- leading terms in 
Hence, the leading term in the quantum component of 
the effective potential only involves the first and second 
derivatives of the potential energy. Since these terms 
already appear in the classical effective potential, eval- 
uating the quantum corrections in this limit does not 
appreciably increase the computational cost of the calcu- 
lation. 

From this point on, the derivation of the DRP for- 
malism with quantum effects is completely analogous to 
the classical DRP approach: the integrand in Eq. (16) 



expresses the statistical weight of the path connecting 
the initial and final configurations, in a time interval t. 
The exponents exp(— S e ff[X]) and exp(— S^[X]) rep- 
resent the classical and quantum contributions to the 
probability of a given path, respectively. In partic- 
ular, the most probable (or dominant) reaction path- 
ways are those which minimize the total effective action 
5[X] = S e ff[X] + S^[X], and are a solution of the 
equation of motion 



1 

2D. 



: d?i = fi(y c //(X) + ^ / (X)) 



with boundary conditions 



X(t) = x /; 
X(0) = x<. 



(22) 



(23) 



The numerical advantage of the DRP approach follows 
from observing that the equation of motion for the dom- 
inant paths conserves an effective energy. In particular, 
Eq. ( 22 ) conserves the quantity 
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4D 
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This property makes it possible to switch from the 
tarae-dependent Newtonian description to the energy- 
dependent Hamilton-Jacobi (HJ) description. To this 
goal, it is convenient to introduce the rescaled atomic 
coordinates 



Vi 



1 



(25) 



where Xi = \fDo/Di is a dimensionless scaling parame- 
ter and .Do is an arbitrary reference diffusion coefficient, 
which has been introduced to ensure the correct dimen- 
sionality of the yi variables. 

Based on this definition, the solutions of Eq 



(22) 



with the appropriate boundary-conditions are obtained 
by minimizing the functional 



Shj[Y] 



i r ' 

DoL d ' iE "> 



y eff [xiyi(i),»,XNyN(i)] -+ 
Veff [xiyi(i),-,XNyN(i)]} 



1/2 



(26) 



where dl = \J Y^=i ^yf • Notice that, since dl oc y/Do, 
yi oc 1/\/Dq and Xi °t V^o^ the arbitrary reference dif- 
fusion coefficient Do cancels out in the effective action 
Shj[x}. 



In the HJ effective action (26) the time variable has 
been replaced by the curvilinear abscissa Z, which has 
the dimension of a length. The crucial point is that in 
molecular systems there is no decoupling of the intrinsic 
length scales. As a result, in order to describe reactions 
as complex as a conformational transition of a peptide 
chain, only about 100 fixed dl steps are usually sufficient 
to reach a convergent discretized representation of the 
integral in Eq. (26). This number should be compared 
with the 10 9 — Kr^MD time steps required to simulate a 
single protein folding transition with mean first-passage 
time in the /is - ms range. 

In the DRP formalism, it is possible to recover the 
information about the real-time evolution of the system. 
In fact the time at which a given configuration of the 
most-probable path is visited is given by the equation: 
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TABLE I: The parameters specifying the two- dim ensional en- 
ergy surface of the toy model defined in Eq. (29). 



In practice, finding the dominant reaction pathway 
amounts to minimizing a discretized version of the ef- 
fective HJ functional: 



n s 



S d H j[Y] = E 1^// + V 'ff 0^)1 A ^.-+i 



*00 = \[^r J Y dl i E eff + Ve f f [xm(i),-,XNy N (i)] 
+v e Q ff [xm(i),-,XNyN(i)]}~ 1/2 ■ (27) 

Notice that also the transition time does not depend on 
the specific choice of the reference diffusion parameter 
Do? as expected. 

The total time is determined by the choice of the effec- 
tive energy parameter E e ff. Its numerical value should 
not be chosen unrealistically large, to avoid introducing 
a bias towards ultra-fast transitions. 



(28) 

where AZ^+i is the Euclidean distance between the slices 

i and i + 1, i.e AZ M+i = \J\Y i+1 - Y^ 2 . 

In the discretized representation of the HJ effective ac- 
tion (28), the width of the distribution of the Euclidean 
distances between consecutive path slices, A/^ +i , should 
not be allowed to increase in an uncontrolled way, in or- 
der to prevent all frames to collapse into the reactant 
or product configurations. As discussed in [23, 24 , the 
most convenient way to achieve this is to introduce a La- 
grange multiplier in the minimization algorithm, which 
holds fixed the ratio between the mean-square deviation 
from the average of the inter-slice distances a 2 of the 
average square inter-slice distance (A/ 2 ). 



IV. ILLUSTRATIVE TEST EXAMPLES 

It is instructive to illustrate the quantum version of the 
DRP approach in simple systems, before tackling realistic 
molecular transitions. 



A. A two-dimensional toy model 

First, we consider the diffusion of a point particle in 
the two-dimensional potential 

3 

U(x,y) = ^Aiexp[-ai(x - x^ 2 - (3i(y - yi) 2 } (29) 

2=1 

The parameters of the potential are given in table [IJ The 
temperature of the heat-bath was set to 300 K, the mass 
of the particle was chosen to be m = 1 u. 

In Fig. [2] we compare the dominant reaction pathway 
obtained in the classical DRP approach (circles) , with the 
one computed keeping into account the quantum correc- 
tion (triangles). In addition we plot the minimum-energy 
path (squares), obtained by minimizing the functional 
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FIG. 2: The dominant reaction pathway in classical Langevin 
dyn amics in the toy model defined by the potential energy 
( 29 ) . The circles denote the classical dominant reaction path- 
way, obtained minimizing the HJ function (26), the squ ares 
represent the minimum-energy path obtained from ( |3Q| ) and 
the triangles the dominant pathway with quantum correc- 
tions, obtained including the quantum component of the ef- 
fective potential in the HJ action. 



which corresponds to the classical dominant path in the 
low-temperature and long transition-time limit [19]. In 
the background we plot the energy map. We observe that 
the quantum corrections on the dominant path are ap- 
preciable. We find that, in this case, it tends to approach 
the minimum-energy path in the transition region. 



B. H2 dissociation on the Cu(llO) surface 

The DRP formalism allows to compute the quantum 
corrections to the real-time dynamics of diffusion- driven 
reactions in non-equilibrium conditions, i.e. for time in- 
tervals much smaller than the thermal relaxation time. 
Quantum effects on reaction kinetics have also been stud- 
ied in the context of instanton-based approaches [5HT2]. 
Such methods are mostly used to evaluate the reaction 
rates and are based on the saddle-point expansion of the 
imaginary-time path integral (i.e. the quantum partition 
function). The corresponding saddle-point paths (i.e. the 
instantons) do not directly relate to physical trajectories, 
hence to the real-time dynamics of the system. However, 
they provide information about the most often visited 
configurations in the transition region at thermal equi- 
librium, and therefore have been used to evaluate the 
change in the free-energy barrier due to quantum effects. 

It is reasonable to expect that the leading quantum 
correction to the free-energy barrier should be qualita- 
tively consistent with the leading quantum correction to 
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FIG. 3: The molecular energy evaluated along the classical 
(circles) and quantum (triangles) dominant reaction path- 
ways, for the H2 dissociation on the Cu(110) surface, at 300 K. 



the energy barrier overcome by the most probable reac- 
tion pathways. Hence, it is interesting to compare the re- 
sults obtained in the DRP and instanton approaches. To 
this end, we consider the H2 dissociation on the Cu(110) 
surface, a reaction which has been investigated in de- 
tail, using instanton methods [9HTT]. In these studies it 
was shown that the quantum corrections lead to an effec- 
tive reduction of the free-energy barrier with respect to 
a classical calculation. In particular, at 300 K the quan- 
tum corrections lower the free-energy barrier by about 
0.1 eV [10]. 

We have used the DRP approach to study the same 
reaction, adopting the same interaction potential (de- 
fined in detail in Ref. [11 J. We have calculated classical 
and quantum dominant reaction paths and used them to 
evaluate the quantum correction to the molecular energy 
barrier overcome along the reaction path. We found that 
at 300K quantum effects lower such a barrier by about 
0.1 eV, which is compatible with the free-energy change 
calculated using the instanton method — see Fig. [3] — . 

A further qualitative insight on the effects of quantum 
corrections can be inferred by comparing the quantum 
dominant reaction pathway and the instanton trajecto- 
ries calculated in [11] . The DRP result is shown in Fig. [4] 
where it is compared to the minimum energy path. We 
see that the dominant reaction path is shorter than the 
MEP, again qualitatively agreeing with the results shown 
in Fig. 3 of Ref. [11]. 



V. CONFORMATIONAL TRANSITION OF A 
PEPTIDE 

We now apply the same method to investigate the role 
of quantum fluctuations in a prototypical bio-molecular 
conformational reaction, namely the C7 eq — » C7 ax tran- 
sition of the alanine dipeptide. 
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FIG. 4: The reaction pathway for the H2 dissociation on the 
Cu(110) surface, projected onto the plane selected by the x 
and y coordinates of one of the two hydrogens. The circles 
denote the minimum- energy path and the triangle represent 
the quantum dominant path. 



Let us begin by showing that the semiclassical approxi- 
mation which underlies the present approach is amenable 
to investigating the conformational dynamics of a pep- 
tide. To this end, we observe that the typical diffusion 
coefficient for an amino acid of mass m — 80 u in water is 
D ~ 1.2 x 10" 3 nm 2 ps _1 , hence 7-6 ps -1 . The 

condition for a semiclassical treatment of the dynamics 
is therefore realized: 7 j3 h c± 0.2. The smallest time scale 
we can reliably describe using the overdamped limit is of 
the order of I/7 ~ 0.2 ps. 

The molecular energy was obtained from the Assisted 
Model Building with Energy Refinement (AMBER99) 
empirical force field [26], without solvent-induced inter- 
actions, at a temperature of 25° C. The initial and final 
configurations of the peptide were obtained by minimiz- 
ing the potential energy. 

For realistic reactions, the global minimization of the 
HJ effective action is in general a challenging task. The 
main difficulties arise from the ruggedness of the ef- 
fective potential and the high dimensionality of molec- 
ular systems. As a result, the most commonly used 
global optimization algorithms — such as e.g. simulated 
annealing — tend to get stuck in secondary minima of 
the action functional. Clearly, in this case, the calcu- 
lated dominant paths would be strongly biased by the 
choice of the initial trial path. 

Our previous tests on molecular reactions have shown 
that the Fast Inertial Relaxation Engine (FIRE) method 
[27] offers a good compromise between performance and 
simplicity [23]. The FIRE algorithm is based on a modi- 
fied dynamics approach and is less prone to remain stuck 
in local minima than other minimization procedures such 
as conjugate gradients or Broyden-Fletcher-Goldfarb- 
Shanno methods. Moreover, we found that the FIRE 
algorithm was more efficient than other global methods 
like for instance simulated annealing. 

The minimization protocol adopted in the present work 




-90 -60 -30 30 60 90 



FIG. 5: The structure of the dominant reaction pathways for 
the the C7 eq —> C7 ax transition of alanine dipeptide, obtained 
in different approaches and projected on the Ramachandran 
plane. In the background is reported the free energy land- 
scape obtained in the classical approach. 



was the following: we generated an initial trial path con- 
sisting of a linear trajectory which connects the initial 
and final points in the Ramachandran plane specified by 
the ip and <\> dihedral angles of the di-peptide pQ . Such a 
path was discretized using 100 equally-displaced frames. 
The path so obtained was initially relaxed the by means 
of a Nudged Elastic Band (NEB) [28] minimization. This 
step is crucial in order to avoid instabilities in the sub- 
sequent DRP minimization algorithm. The NEB path 
was then used as a starting point for the minimization of 
the classical DRP action, followed by the minimization of 
the complete quantum DRP action. The effective energy 
parameter E e ff was chosen to be 10% larger than the 
maximum value of along the NEB path. This 

condition ensures a long transition time, and avoids that, 
during the minimization, the HJ effective action becomes 
complex. 

In Fig. [3] we plot the dominant paths obtained in the 
DRP approach with and without quantum corrections, 
together with the minimum-energy path. In the back- 
ground we show the free-energy as a function of the di- 
hedral angles, evaluated by means of an all-atom classical 
meta-dynamics simulation [29] [30], using the same force 
field. In Fig. [6] we present the evolution of the molecular 
energy along the reaction path, while in Fig. [7] we re- 
port the evolution of the distance between the H18 and 
the 06 atoms, which are involved in a hydrogen bond. 
The entire set of calculations required in total about 700 
hours on 2.2 GHz processors. 

Some comments on these results are in order. First 
of all, we observe that the classical dominant pathway 
crosses the barrier in a region in which the molecular po- 



tential energy is about twice as large than at the saddle 
point, which is visited by the minimum-energy path, see 
Fig. [6j We emphasize the fact that the dominant reac- 
tion pathways are expected to describe genuinely non- 
equilibrium transitions. In general, such paths do not 
need to cross the barrier precisely at the saddle-point. 
On the other hand, it is important to check that such 
a large effect is not an artifact of the calculation. For 
example, problems may emerge if the effective energy 
parameter E e ff was chosen very large. In this case, the 
total transition time would be very small and the cal- 
culation would lead information about the dynamics of 
ultra-fast transitions. In addition, problems may emerge 
if the path space was not sufficiently explored during the 
minimization procedure. In order to check the numerical 
reliability of our results, we have computed the classical 
dominant pathway, starting from a path which crosses 
the barrier at the saddle point, with an effective energy 
only 1% larger than the maximum value of \V e f f [x] | along 
the initial path. After the minimization of the HJ action, 
we recovered the same result for the classical dominant 
path shown in Fig. |3j and a very similar transition time. 
This result makes us confident that the dominant paths 
are independent on the choice of the initial trial path and 
are not appreciably dependent on the specific choice of 
the effective energy. 

A second important result of our calculation is that 
the quantum effects on the structure of the dominant 
reaction pathways are clearly visible. Even though all 
the three paths are qualitatively similar, the energy dif- 
ference which is overcome by the most probable reaction 
pathway in the presence of quantum fluctuations is about 
50% smaller than the one in the classical case (see Fig. [6]). 
Such a large difference arises because the molecular en- 
ergy surface in the transition region is quite steep and 
the quantum fluctuations of hydrogen atoms are quite 
large. Note that the (classical) free energy differences in 
the same region are much smaller, due to a relatively high 
entropic contribution, associated e.g. to the rotation of 
the other dihedral angles. In addition, the distance be- 
tween the atoms involved in the hydrogen bond is always 
about 0.2 A larger in the quantum than in the classical 
dominant path (see Fig. [7|. This suggests that the en- 
ergy difference overcome by the classical dominant path 
is larger, because the O and H atoms get closer, hence 
increasing their van der Waals repulsion. 

In Fig. [8] we compare the time evolution of the system 
in the classical and quantum calculations, i.e. we plot 
the time at which each value of the reaction coordinate 
is visited along the transition. First of all, we note that 
in both calculations the most probable transition lasts 
about 15 ps, i.e. a time much longer than the 0.2 ps time 
scale below which the overdamped approximation is no 
longer appropriate. It is also interesting to compare the 
velocity of the classical (circles) and quantum (triangles) 
dominant transition, which is represented by the slope 
of the curves l(t). We see that while the velocity along 
the classical path is essentially constant throughout the 
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FIG. 6: The evolution of the molecular energy along the re- 
action coordinate / of the C7 eq — ► C7 ax transition of alanine 
di-peptide. The circles denote the classical dominant reac- 
tion pathway, obtained minimizing the H J function ( 26 ) , the 



squares represent the minimum-energy path obtained from 
(30) and the triangles the dominant pathway with quantum 
corrections, obtained including the quantum component of 
the effective potential in the HJ action. 
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FIG. 7: The evolution of the distance between the H18 
and the 06 atoms, along the reaction coordinate / of the 



C7 e 



CI ax transition of alanine di-peptide. The circles de- 



note the classical domin ant reaction pathway, obtained min- 
imizing the H J function ( 26 ) and the triangles the dominant 



pathway with quantum corrections, obtained including the 
quantum component of the effective potential in the HJ ac- 
tion. 



entire reaction, the quantum dominant path accelerates 
after about 3 ps after about 4 ps, i.e. in the region of 
high force before and after the saddle. 



VI. CONCLUSIONS 

In this work we have introduced a formalism which al- 
lows to investigate at the semi-classical level the role of 
quantum fluctuations of atomic nuclei in the real-time dy- 
namics of non-equilibrium molecular transitions. Unlike 
other method which are more suited for rate calculations 
[SJ [TTJ [12] or exploring the short-time dynamics inside 
a thermodynamical state [7], the present DRP approach 
is particularly efficient in investigating the real time dy- 
namics as the system is crossing the free-energy barrier. 
The computational efficiency of the method makes it pos- 
sible to study reactions involving large molecules, such as 
e.g. peptide chains. 

From an analysis of the quantum corrections to the 
Langevin equation we have shown that the quantum 
stochastic dynamics in the vicinity of the saddles with 
negative Hessian trace is qualitatively similar to a clas- 
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FIG. 8: The time evolution of the dipeptide, in the C7 eq 
CI ax transition of alanine di-peptide. The circles denote the 
classical domi nan t reaction pathway, obtained minimizing the 
HJ function (26), the triangles the dominant pathway with 



quantum corrections, obtained including the quantum com- 
ponent of the effective potential in the HJ action. In the in- 
serts we show the position of some frames along the reaction 
path. 



sical stochastic dynamics at a lower temperature. Con- 
versely, in the vicinity of saddles with positive Hessian 
trace, or in the potential wells, quantum fluctuations can 
be interpreted as effectively raising the temperature. 

We have shown that, in the test case of H 2 dissoci- 
ation on copper, the quantum corrections obtained in 
the DRP formalism away from equilibrium qualitatively 
agree with those obtained in equilibrium conditions, us- 
ing the instanton method. 

We have applied the DRP formalism to the study of 
the C7 eq — » C7 ax transition of alanine dipeptide which 
represents a prototypical example of biomolecular tran- 
sition involving a hydrogen bond. We have found that in 
this reaction the inclusion of quantum fluctuations can 
significantly modify the reaction path with respect to a 
classical calculation. 

We conclude this work by discussing possible limita- 
tions of the present approach. In general, we expect the 
DRP method to become inefficient in the following sce- 



narios: 



• For each different boundary condition (23) there 
exists a large number of local minima of the HJ 
functional, all with comparable statistical weight, 
exp(-Sjyj). 

• The reaction mechanism depends very strongly on 
the initial configuration and the reactant space 
is large. In this case, a very large number of 
reaction pathways would be needed in order to 
fully characterize the transition. By contrast, any 
method which provides only a relatively small num- 
ber of them would carry insufficient information 
(note that this limitation applies also to MD simu- 



lations). In this case, one must rely on a description 
based on the projection on a small set of reaction 
coordinates. 

• The fluctuations around each of the different dom- 
inant paths are very large and the regions visited 
by the fluctuations associated to different dominant 
paths significantly overlap. In this case, the very 
notion of dominant pathway looses its significance. 
On the other hand, if such fluctuations are rela- 
tively small, their contribution can be systemati- 
cally included through a perturbative expansion in 
the thermal energy /c#T, using the method recently 
developed in Ref. [TB] . 

In order to assess how such potential limitations affect 
the applicability of the DRP method to realistic molecu- 
lar transitions, several comparative tests were performed, 
based on the comparison against the result obtained 
by MD simulations [2TJ [22] . These studies have shown 
that the DRP approach yields the correct description of 
the non-equilibrium dynamics of complex macromolec- 
ular transitions, such as protein folding. On the other 
hand, we emphasize that the semi-classical approach pre- 
sented here works in conditions in which the quantum 
effects provide at most small corrections to thermally ac- 
tivated pathways, hence in the presence of dissipative dy- 
namics. It is not applicable to investigating the non- 
dissipative tunneling and in general the dynamics in the 
deeply quantum regime. 
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Appendix A: Multi-dimensional Generalization of 
the Quantum Smoluchowski Equation 

In this appendix, we provide the generalization of the 
QSE Q to a physical system consisting of 3 N degrees of 
freedom (e.g. the atomic coordinates of a molecule). We 
begin by observing that the naive substitution ^ — > V 
in Eq. ^ does not represent the correct generaliza- 
tion, as it does not lead to the correct equilibrium distri- 
bution P eq (x), describing the thermodynamical limit in 
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the semi-classical regime. In order to obtain the correct 
multi-dimensional generalization of the QSE ([7) one can 
use the same path integral approach adopted in]4 . Since 
such a procedure is quite lengthy and technically rather 
involved, here we present an alternative, albeit slightly 
less rigorous, derivation leading to the same result. 

The conservation of the number of particles implies a 
continuity equation for the probability, i.e. 

N 

3 t P(X,*)=V.J(X,*) = ^Vi • Ji(X,t), (Al) 



where J(X,t) = (jx(X, f), j 2 (X, t), . . . , ^(X, £)) is the 
probability current. 

Without loss of generality, each single-particle compo- 
nent of the current ji(X, t) can be defined in terms of a 
set of in-so-far unspecified functions £i(X), and QOQ'- 

f{X,t) = Dl (/3V i C/(X)+el(X))p(X,t) 

+ V*[£>* (l + £(X))P(X,t)], (A2) 

where D l = l/(rriif3j) are the classical diffusion con- 
stants. Such a definition assures that in the limit 
£|(X),£2pO ~* 0? one recovers the classical Smolu- 
chowski equation. The functions 

D\(X) = Di (/3V i f/(X)+el(X)) (A3) 

D\{X) = Dl (1 + £(X)) (A4) 

are the multi-dimensional generalization of the Moyal co- 
efficients discussed in [4]. 

The unknown functions £{(X) and ClPO c an b e de- 
termined by requiring that the continuity Eq. (Al) must 



Hence, 



yield the correct thermodynamics, i.e. that its station- 
ary solution coincides with the well-known semi-classical 
expression of the Boltzmann's weight, 

P eq (X) = exp(-/?C/(X)) (l-i 1 (X)+L 2 (X)), (A5) 

with 

N 



MX) 
MX) 



(3j2\ k V 2 fe P(X) 



(A6) 



k = l 

N 



P 2 



J2^k |V,/7(X)| 2 (A7) 



k=i 



Using Eq. (A2) and Eq. (A5) and imposing the 



condition of vanishing current at thermal equilibrium, 
limt-^oo j(X, t) = 0, up to leading order in the quantum 
expansion parameters, we obtain 

£l(X) + V^(X) + V%(X) 

-V'Li(X) - ^(X)/?W(X) = (A8) 

The consistency with one-dimensional result ^ implies 

6(X) = Li(X). (A9) 



£(X) = L 2 (X)/3W(X) - V%(X). (A10) 



The QSE obtained from Eq.s (jATJ, (|A2f , (|A9) and ( |A1Q| ) 
is equivalent to Eq. (J8|, to leading order in the quantum 
expansion parameters However, the form ([8) is usu- 
ally preferred, as it guarantees the consistency with the 
second law of thermodynamics [5]. 



Appendix B: Path Integral Representation of 
Quantum Langevin Dynamics 

In this appendix, we show that the solution of the QSE 
(|8| can be sampled by integrating an associated quantum 
Langevin equation, with a multiplicative noise. We also 
construct the path integral representation (16), which is 



used to derive the quantum extension of the DRP for- 
malism. 

Let us begin by considering a generic Langevin equa- 
tion with multiplicative noise, in the form 

k i = f i (X)+g(X)r} i (t), (i = l,...,N) (Bl) 

where X{ denotes the coordinates of the i— th particle, 
ffi(x) is a 3-dimensioanl stochastic force of unit variance. 
Such a stochastic differential equation generates a prob- 
ability distribution which obeys the generalized Smolu- 
chowski equation — see e.g. discussion in [25] — 



d_ 

at 



P(X,t) = ^V, [(-^(X)-^(X)V^(X)) P(X,t) 

(B2) 



\v t ( 5 2 (X) P(X,t)) 



The real parameter < a < 1 specifies the s toch astic 
calculus adopted to define the differential Eq. (Bl). In 
particular, a = (a = 1/2) corresponds to the so-called 
Ito (Stratonovich) calculus. 

We now want to derive a Langevin equation in the form 



(Bl) which generates a probability density obeying the 
QSE (J8|. To this end, it is important to emphasize that 
QSE is an ordinary partial differential equation, hence it 
is uniquely defined in the standard (e.g. Riemann) cal- 
culus. Hence, for every choice of stochastic calculus a 
there is in general a different Langevin equations, asso- 
ciated to the same physical Smoluchowski Eq. This 
can be obtained by finding the functions ^(X) and the 
vector fields /i(X) such that (B2) Smoluchowski equation 
coincides with the QSE ([8]), to order A accuracy. Such a 
request leads to 



5i (X) = v/2A ^l + -Li(X)J (B3) 
fi(X) = -DiPViU(X) + DiP$i(X), (B4) 
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where Qi(X) is a "quantum force" whose definition de- 
pends on the calculus adopted and reads 

Qi(X) = ^V,L 2 (X) 

- L 1 (X)V i f/(X)--V i L 1 (X). (B5) 

Also the path integral representation of the solution 
of the quantum Smoluchowski Eq. depends on the 
choice of the stochastic calculus and reads [25] 



Note that the path integral now contains the standard 
Wiener measure 



N t N 



vx = hm nn 

AT* -yon -L J- J- J- 



dxi(l) 



jv^oo 1111 UDiirAt) 3 / 2 ' 

i=i %=i v / 



(B9) 



P(X,t|X 



PX exp 



i=l 



1 



2 5l 2 (X) 



ii - /i(X) + a ff (X)V i < ?i (X) + aVi • / 4 (X) 



}■ 

(B6) 



The first exponent can be taken out of the path inte- 
gral since it does not affect the statistical weight of the 
diffusive paths. To see this, we introduce a scalar func- 
tion W{x) which is defined as the formal solution of the 



where <fc(X) and /»(X) are given by Eq.s (|B3j) and (|B4j), partial differential equation 
and the modified Wiener measure VX depends on the 
configuration and reads 



N t N 

VX = lim WW 



dxi(l) 



\t =1 mnAt (1 + A^V?C/[X(Z)])] 



3/2 ' 



V*W(X) = (-V^(X) + (3Q l - aW l L 1 (X)^j (BIO) 
(B7) With such a definition, the first exponent in Eq. ( |B8| ) is 
written as an exact differential form, 



where At = t/N t . Plugging p3|)-p4|) into Eq. p6 ) 
and expanding to leading order in the we obtain, after 
some tedious but rather straightforward calculations 



P(X,*|Xi 



VX exp 



-W-aV^lX))]} 

A/3 2 



exp 



2~ 



(/3Vi£/(: 



4A 



4 vl Vi£/(X)| 2 - |v?£/(X)) + apDiVi ■ &(X) 



-^/3 2 |V i f/(X)| 2 L 1 (X) + a/3AViii(X) • V;C/(X) 



= e -f(^(X / )-lV(X l )) 

= ^(X/.Xi) (Bll) 



L . which depends only on the end-points and not on the 
J path. If we now specialize on the Stratonivich calculus 
(B8)we obtain Eq. (fiB. 
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